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A hierarchical Bayesian setting for an inverse 
problem in linear parabolic PDEs with noisy 
boundary conditions 

Fabrizio Ruggeri * * * § , Zaid Sawlan t, Marco Scavino * and Raul Tempone § 

Abstract. In this work we develop a Bayesian setting to infer unknown parameters 
in initial-boundary value problems related to linear parabolic partial differential 
equations. We realistically assume that the boundary data are noisy, for a given 
prescribed initial condition. We show how to derive the joint likelihood function 
for the forward problem, given some measurements of the solution field subject to 
Gaussian noise. Given Gaussian priors for the time-dependent Dirichlet boundary 
values, we analytically marginalize the joint likelihood using the linearity of the 
equation. Our hierarchical Bayesian approach is fully implemented in an example 
that involves the heat equation. In this example, the thermal diffusivity is the 
unknown parameter. We assume that the thermal diffusivity parameter can be 
modeled a priori through a lognormal random variable or by means of a space- 
dependent stationary lognormal random field. Synthetic data are used to test the 
inference. We exploit the behavior of the non-normalized log posterior distribution 
of the thermal diffusivity. Then, we use the Laplace method to obtain an approx¬ 
imated Gaussian posterior and therefore avoid costly Markov Chain Monte Carlo 
computations. Expected information gains and predictive posterior densities for 
observable quantities are numerically estimated using Laplace approximation for 
different experimental setups. 

Keywords: Linear Parabolic PDEs, Noisy Boundary Parameters, Bayesian Infer¬ 
ence, Heat Equation, Thermal Diffusivity. 


1 Introduction 

Parabolic partial differential equations model various important physical phenomena 
such as diffusion and heat transport. The solution of such equations propagates forward 
in time from an initial condition given boundary conditions and equation coefficients. 
In applications, some equation coefficients can be unknown quantities that need to be 
estimated. In addition, exact initial and boundary conditions might not be known. One 
possible approach to estimate these unknowns is to solve the inverse problem given some 
information about the solution. Classical inversion methods for parabolic partial differ- 
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ential equations are introduced by Samarskii and Vabishchevich (2007) in the Chapter 
8 of their book. 

In this work, we consider a Bayesian inversion problem to determine the coeffi¬ 
cients of linear parabolic partial differential equations, under the assumption that noisy 
measurements are available in the interior of a domain of interest and for the unknown 
boundary conditions. A main novelty of our approach to solve the inverse problem relies 
on the assumption that the boundary parameters are unknown and modeled by means 
of adequate probability distributions. Subsequently, the contribution of the boundary 
parameters is marginalized out from the joint law with the unknown equation coeffi¬ 
cients we want to infer, allowing the characterization of their posterior distribution. 
There are many advantages to a Bayesian approach. For example, it provides a solution 
along with a comprehensive measure of uncertainty given by the posterior distribu¬ 
tion. Moreover, the prior available information can be easily incorporated in terms 


of elicited prior distributions (Ghosh et al. 2006). An important issue in this work is 
that Bayesian inversion is posed as a hierarchical process. Boundary conditions can 


therefore be treated as unknown parameters (Kaipio and Fox| [2011). Since we are only 
interested in estimating the equation coefficients, we eliminate those extra parameters 
by marginalization. 

Bayesian inversion techniques for the heat equation have been discussed and imple¬ 
mented in some previous works. Kaipio and Fox (2011) provided a general Bayesian 


framework for inverse problems in heat transfer, classified according to the dominant 
mode in heat transfer. The authors address many issues regarding forward problems 
and their statistical analysis. The prior modeling is extensively discussed, as well as 
how to deal, in particular, with different sources of uncertainties. Heat flux reconstruc¬ 


tion problems have been studied by Wang and Zabaras (2005 2004). When referring to 


the problem of parameter estimation in inverse heat conduction problems, Wang and 
Zabaras showed how to infer the thermal conductivity using a hierarchical Bayesian 
framework, on the basis of temperature readings within a conducting solid, assuming 
that the heat flux on the boundary and the heat source are known. They also explored 
the high dimensional posterior state space by means of Markov Chain Monte Carlo 


simulation. Lanzarone et al. (2014) estimated the thermal conductivity of a polymer, 


transforming the heat equation into a stochastic differential equation and considering the 
Euler-Maruyama approximation to get the likelihood, introducing latent observations 
in space and then using a relatively cumbersome Markov Chain Monte Carlo method. 


Fudym et al. (2008) and Massard et al. (2010) addressed the estimation problem of 


the thermal diffusivity, as in the present work, which is a parameter that describes 
thermophysical property of materials. In their works a large number of temperature 
measurements is made by an infrared camera, with fine spatial resolution and high 
frequency. They solved one and two-dimensional forward problems for transient heat 
conduction, with spatially varying thermal conductivity and volumetric heat capacity, 
by finite differences, according to a nodal strategy. The parameter vector at each node 
is then estimated either by minimizing an a posteriori objective function when prior 
Gaussian distributions are assumed for the parameters, or by means of Markov Chain 
Monte Carlo methods for different prior distributions. 
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This work is organized as follows. In Section [2j we introduce the statistical setting 
and we derive the explicit form of the joint law of the unknown equation coefficients and 
the boundary parameters. In Section [3j we use a finite element scheme in order to write 
the solution of the forward problem as a linear function of the boundary conditions. 
We demonstrate in Section [I] that, under certain conditions, an exact marginalization 
can be carried out, yielding a closed formula for the marginal likelihood of the equa¬ 
tion coefficients. In Section [5j we apply our marginalization technique to estimating 
thermal diffusivity in the one-dimensional heat equation in two cases given temperature 
simulated data. Numerical results are obtained for the non-normalized log posterior 
distribution of the thermal diffusivity. We model prior knowledge about the thermal 
diffusivity first as a lognormal random variable and then using a lognormal random 


field with a squared exponential (SE) covariance function (Rasmussen and Williams 


2006). In the first case, we use the Laplace method to provide an approximated Gaus¬ 


sian posterior distribution for the thermal diffusivity. Such method is then applied to 
obtain fast estimations of the information gain and the expected information gain under 
three experimental setups, and the predictive posterior mean of the temperature is also 
derived using the inferred thermal diffusivity. In the second case, where the thermal 
diffusivity is allowed to depend on the spatial variable, the Laplace approximation is 
used to obtain the posterior distribution of the hyperparameters that characterize the 
prior distribution for the thermal diffusivity. 


2 Statistical setting and preliminary results 

In this section we introduce the statistical model associated to the forward initial¬ 
boundary value problems for linear parabolic partial differential equations. We then 
derive, under mild assumptions, the exact expression for the joint likelihood function of 
the unknown parameters in the parabolic equation and the unknown boundary param¬ 
eters. 


Consider the deterministic one-dimensional parabolic initial-boundary value prob¬ 
lem: 


( d t T + L e T = 0, 

I T(x l , t) = T L (t), 
I T(x R ,t) = T R (t ), 
0) = g{x ), 


x £ (xl,x r ) : 0 < t ^ tw < oo 
t £ [0, f at] 
t £ [0, t N ] 
x £ ( x L ,x R ), 


( 1 ) 


where Lg is a linear second-order partial differential operator that takes the form 
L e T = -d x {a{x)d x T) + b{x)d x T + c{x)T , 

0(x) = (a(x), b(x), c(x)) tr , and the partial differential operator d t + Lg is parabolic, 
because (Evans (1998), [p.372]) there exists e such that a(x) ^ e > 0 for all x £ {xl,x r ). 
We also assume that 


PI a, b and c are bounded functions. 







4 


Bayesian Inference for Linear Parabolic PDEs 


P2 TljTr and g are square integrable functions. 

P3 The initial condition, g , is consistent with the boundary functions, namely g{x R ) = 
Tl( 0) and g(x R ) = T R ( 0). 


Then, under the assumptions P1-P3, there exists a unique weak solution of 0 0 vansj 
(1998]), [pp.375-377]). 

Our main objective is to provide a Bayesian solution to an inverse problem for 0, 
where we assume that 


i Q is unknown, while the initial condition g in the initial-boundary value problem 
is known; 

ii Q is allowed to vary with the spatial variable x. 


Remark 2.1 In our Bayesian approach, we will assume later that the coefficient a(x) 
is a lognormal random variable or lognormal random field. Therefore, a(x) will not be 
bounded as assumed in PI. However, it can be proved that there exists a unique solution 
of the stochastic parabolic initial-boundary value problem in the space L 2 {fl, H 1 ). Such 
proof can be found in Charrier (2012) for elliptic boundary value problems but it can be 
also extended to parabolic initial-boundary value problems. 


Given noisy readings of the function T(x,t) at the I + 1 spatial locations, including 
the boundaries, Xl = Xo,Xi,x 2 , ■ ■ ■ ,xi-\,xi = x R , at each of the N times t\,t 2 , ■ ■ ■, ijv, 
we want to infer 6 using a Bayesian approach. To determine the posterior distribution 
for 0, we need first to obtain the likelihood function of 0. The remainder of this 
section derives the joint likelihood function of 0 and the boundary parameters. Let 
us introduce some convenient notation and assumptions: let Y n := (Yo,n; ■ • • i Yi <n ) tr 
denote the vector of observed readings at time t n , and assume a statistical model with 
an additive Gaussian experimental noise e n ; that is: 


Y n 

(/ + l)xl 


T L (t n ) 

T(xi,t n ) 


T R (t n ) 




n ■> 


( 2 ) 


where e n *~ d ' Af(0 I+1 , a 2 I/+i) for some measurement error variance a 2 > 0. The 
covariance matrix of e n is assumed equal to er 2 I/+i for simplicity, a general covariance 
matrix E £tj could be used as well provided that the boundary measurement errors are 
independent from the interior measurement errors. 

Also denote by := (Yi t „,..., Y/_i in ) tr the vector of observed data at the 

(J-1)X1 

interior locations x\,x 2 ,- ■ ■, aq_i and let Y^ := (Y L , n , Y RyU ) tr be the vector of observed 
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data at the boundary locations Xq,Xj at time t n . 

The density of Y* is derived as it follows. First consider the time local problem, defined 
between consecutive measurement times, i.e. 

{ d t T + L g T = 0, x £ (x L ,x R ), t „-1 <t^t n , 

T(x L ,t) =T L (t), t€[t n -i,t n \, 

T{xr, t) = T R (t), t G [t n -i,t n \, 

T{x , t n —\) T(x, t n — i), x G ( x R , Xr ), 

whose exact solution, denoted by T(-,t n ), depends only on 9, T(-,t n - 1 ) and the bound¬ 
ary values {Ti,(t), Tj{(t ) tn y Finally, use the form of the statistical model |2| to 
obtain the form of the density of Y*. 

Lemma 2.2 Given the model |i|), the probability density function of Y* is expressed 
as 

P(Y*|0,?(-,*„_!), {T L (t),T R (t)} Kltn _ i>tn) ) = lIRtJI^) , (4) 

where := {T(xi, t n ) — Yi, ra , ..., T(xj_i, t n ) — Y/_i ira ) tr denotes the data residual 

{I 1) X 1 

vector at time t = t n . 


For illustration purposes and without loss of generality, assume now that the Dirich- 
let boundary condition functions, Tr(-) and Tr(-), are well approximated by piecewise 
linear continuous functions in the time partition {t n } n=1 N . 

In this way, only 2N parameters, say T L (t n ) = T Li71 , T R (t n ) = T R>n , n = 1, 2 
suffice to determine the boundary conditions that are, in principle, infinite dimensional 
parameters. Let LBC n denote the time nodes that determine the local boundary con¬ 
ditions {T£,, n _i, Ti )T j, Tr iTI —i, Tr h } . Other interpolation schemes may be used as well. 


Remark 2.3 Given the discretized Dirichlet boundary conditions introduced above, we 
can say that T(-,t n ) depends only on 9, T(-,t n - 1 ) and the boundary parameters LBC n . 
Similarly, T(-,t n - 1 ) depends on 9 ,T{-,t n - 2 ) and the boundary parameters LBC n -\. 
From this recursion, we can obtain 

p(Yi|0, 5 ,{LR^}. =1 ; J = p(Yi|0,T(.,t„_ 1 ),LRC'„). (5) 

Since the initial condition, g, is assumed to be known, it will be omitted in the rest of 
the paper. 


Lemma 2.4 Given the model Q) and Lemma 2.2. the joint likelihood function of Q and 
the boundary parameters {LBC n } n _ 1 N is given by 


1 /I \ 

p(Yi,..., Y n \e, {LBC n }„, . „) = n exp ( - ^ ||R,„ ft J 

x 7^ exp (~ h {Tl •" - n ’“ )2 ) x 7^ “ p (“ ^ (Tr “ ■ Yx - nf ) (6) 
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Proof Observe that Y*, Y® are conditionally independent given 9, T(-, t n _i), LBC n . 
Thus, we have 

p(Y n \9,T(;t n ^),LBC n ) = p{Yl,Yl\d,f{-,t n ^),LBC n ), 

= p(Y*|0, !?(.,*„_!), LBC7 n ) x p(Y*\0,f{;t n ^),LBC n ), 
= p{Yl\e,f{-,t n ^),LBC n ) x p(Y*\LBC n ) 

since Y„ | LBC n does not depend on either 9 nor T(-, 

The joint likelihood function can then be written as 

p(Yi, • • •, Y n |0, {LBC n } n=1 n ) 

= p( Y N |0, {LBC n } n=1 _ N , Y N _i,.. •, Y 1 ) x p(Y u ..., Y^O, {LBC n } n=1 ^ N _ l ) 
(since Yn|0, {LBC„} n= i N is independent from Yi,..., Yn-i) 

= P(Y N \0, {LBC n } n=1: n ) x p(Y 1; ..., Yn.,10, {LBC n } n=lt n _ x ) 

(using equation ([5])) 

= p(Y N |0,f (-^jv-r), LBC n ) x p( Y 1( ..., Y N _ X | 9. {LBC n } n=1> 

(iterating the previous arguments) 

N N 

= '[[p(Y a \0,f(-,t n - 1 ),LBC n ) = '[[p(Y I a \0,T(-,t n - 1 ),LBC n ) x p(Y*\LBC n ) 

n =1 n =1 

and finally, using Q, the expression © is obtained. j§ 

Remark 2.5 A generalization of the joint likelihood can be obtained given serial corre¬ 
lations, that is, {e„}^L 1 are time correlated. 


Using the notation T l = (T L>1 ,..., T L , N ) tr , T R = (T R p,... ,T R>N ) tr ,Y L = (Y L<1 ,..., 
and Yr = (Yr i, ..., Y R ^) tr , the joint likelihood function ([b]) can be written as 


/ 1 N 

p( Yr, • • •, Y n |0, T l ,T r ) = (V^a)~ N ^ exp ( ^ ||R, 




x exp 


1 

2^ 2 L 


iTi-Yill^ + llTfl-Yflll^ 


( 7 ) 


which is a suitable expression to derive later the marginal likelihood of 6. The prior 
distributions for 9 and the boundary parameters Tl,T# can be specified in different 
ways. Generally speaking, the prior distribution for 9 should be proposed according to 
the physical properties described by the unknown parameters and taking into account 
available prior knowledge about the observed phenomena. As an example, it is known 


Y l , 
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that the thermal conductivity of polymethyl methacrylate (plexiglas) is in the range 
0.167-0.25 W/(mK). A uniform prior distribution for the thermal conductivity could be 
chosen if there is no preference among the values of the interval. 

As per the boundary parameters, we may assume independent Gaussian prior distribu¬ 
tions Tl j71 ~ A r(pL,n,o'p), Tr^ ~ A crp), n = 1 that are centered at a 

least square spline fit, say p L = {p L ,i, ■ ■ ■, RL,N) tr , Rr = (rra, ■ ■ ■, RR,N) tr , of the 


JVxl 


JVx 1 


observed data, and "Yr, respectively, and for some prior variance > 0. 

We claim that the data residual vector can be written as a linear function of TA 
and Tfl. The next section is devoted to the proof of this basic result. This proof allows 
for the exact marginalization of the contribution of the nuisance boundary parameters 
from the joint likelihood function (J7|. 

3 Numerical Approximation 

In this section, our goal is to approximate the residual vector R tn as a linear function 
of the boundary conditions. After reformulating the main problem ([!]), as described in 
the next lemma, we will introduce its weak form and finite element method will be then 
applied to provide a numerical approximation of the solution of the weak problem. 

Lemma 3.1 The solution of (JT]) can be written in the form 


T(x,t ) =T L (t )——— +T R (t )——— +u(x,t) 


( 8 ) 


Xr - X L Xr- X L 


where u solves a new initial-boundary value problem with homogeneous Dirichlet bound¬ 
ary conditions: 


d t u + Lgu = f(x, t), x £ (xl, xr),0 < t < fjv < oo 
u(x L ,t) = 0, t€[0,t N ] 

u(x R ,t) = 0, te[0,tjy] 

u(x,0) = g°(x), x<=(x l ,xr) 


(9) 


and 



We now introduce the weak formulation of problem Q. 
Find u(t) GV = Hq(xl,xr), t £ (0,tjv) such that: 


Ixl d tu{t)vdx + B ( u(t),v) = f** f(f)vdx,Mv £ V, t £ (0, £jv), 
w(0) = g°, x £ (x L ,x R ), 


(10) 
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where B(u,v) = f^[a(x)d x ud x v + b(x)d x uv + c(x)uv] dx and Hq(x r , x R ) is the closure 
of the space C^(xl, x r ) of continuously diff erent iable functions with compact support 
on (xl,xr) with respect to the Jf^norm (Johnson (1987), [p.149]). 


Given a mesh xl = Xq < ... < kAx = x k < ... < IAx = Xi = x R of the spatial 
domain (x R ,x R ), we apply the finite element method with piecewise linear functions 
(hat functions) {<fk}k=\ to approximate the weak solution u(x,t ) of as linear com¬ 
binations of the basis functions: 


i -1 

u(x , t ) « UAx (*,t)=5> ( t)(j)k (x) , 0 < t ^ t N , 
k= 1 


and we get 
i -1 


i-i 


/ xr M " rXR 

(t> k (t)jdx+^u k {t)B (cfrkAj) = / f(t)<j>jdx,j = l,...I-l,t£(0,t N ). 

K = 1 L k= l ^ XL 

( 11 ) 

Let u(f) = ..., , then we can write the linear system of ODEs 

in a matrix form as: 

Md t n(t) + Sgu(t) = f(f), 0<t.^tN, (12) 

where M is the mass matrix, Sg is the stiffness matrix and f(t) is the load vector. 

We consider now a uniform time discretization of (OAn) such that to = 0 < ••• < 
nAt = t n < ... < NAt = fjv and denote u(f n ) = u n and f(t n ) = f„. We apply the 
backward Euler method on equation (12) to obtain the fully discrete analogue of (10) 
that takes the form 


(M + AtS e ) u„ +1 = M u„ + Atf n+1 , n = 0, ..., N - 1, 
Mu 0 = g°, 


(13) 


where g° = (/** g°(fidx,..., f** . 

Theorem 3.2 The approximation of the weak solution of ([9]) can be written as a linear 
function of the initial and boundary conditions: 

u n = A n (0)\l 0 + A R ^ n {G)L R + A R ^ n (d)P R , (14) 

where u„ = (uq„,..., 'U/- 1 , n ) tT \ T L = (T Lj1 , ..., T LiN ) tr , T R = (T R}1 ,..., T RiN ) tr , 
and the matrices A n (d), A R ^ n (9) and A R , n (Q) are explicitly constructed in the proof. 


Proof See Appendix A <©■ 
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Theorem 3.3 The approximation of the weak solution of |l]) can be written as a linear 
function of the initial and boundary conditions: 

T„ = B"T 0 + A L , n {0)T L + A R , n (0)T R . (15) 

where T n is defined similarly to u ra and the matrices B, A R , n {0) and A R ^ n {6) are ex¬ 
plicitly constructed in the proof. 

Proof See Appendix B (Hi- 

Corollary 3.4 The data residual vector R tji = ( T(x\,tn ) — Y . n , T(x/_i,i n ) — 
Yj—i t n) tr is approximated by 

R= (B"T 0 - Y*) + A L , n (0 )T l + A Rtn (0)T r . (16) 

4 The marginal likelihood of 6 

Given the observations Yi,...,Yn, we showed, in the Section [5] how to obtain the 
joint likelihood function of 0 and the boundary parameters Tl, T fl . 

In the present section, we derive a convenient expression for the marginal likelihood 
of Q, under the assumption that the prior distributions for Tl and T/? are independent 
Gaussian: 

T l ~ Af(/x L ,(7pIjv), 

Tr ~ A/"(Mr> crpljv), 

where p, L and fi R are least square spline fits of the observed data, Yl and Yr, respec¬ 
tively. Then, using 0. the marginal likelihood of 6 is given by: 



X6Xp (~2^ (Ti " - Y l ) - ^2 (T* - YrHTr - Y R )j 

XeXP ( _ 2^ (Ti “ V L ) tr (TL - Ml) - 2^2 ( T r - p R ) tr (T R - p R )j dT L dT R . (18) 


To provide the exact expression of the marginal likelihood of the linear parabolic 
equation coefficients d, which has been implemented in the computational examples 
presented in Section [5j it is convenient to introduce the following notation: 

N N 

Al = y~] A L ^ n (d) tr A L}n (G), A r = A R ^ n {6) tr A R ^ n (0 ), 

NxN , NxN , 

71=1 n= 1 
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N 


N 


A 2jL = ]T A L , n (0) tr (Yi - B”Tq) , A 2 ,* = £ A R , n {0) tr (Yl - B' l T 0 ), 


Nx 1 


n= 1 


Nxl 


n—1 


Alb. = Y A Lyn (0) tr A R , n (G) , D a 2 = diag (\ ( ... A-\ £) , = diag f -4 > • • • > 4 ) • 

iVxiV ^ NxN w CT-J Nx ? N \ ft* I 


Theorem 4.1 The marginal likelihood, of 0 is given by: 

p(Y ll ...,Y N \9) = (V^a)- N ^+ 1 \V^a p )- 2N (2n) N / 2 \A Q \ 1 / 2 (2n) N / 2 \A 1 \ 1 / 2 

N 


1 


1 


x exp } 2^2 [vYlJ-L + /*«*>*] “ 2^2 


1 


YffY l + Y£Y h + ^(Y* - B"T 0 ) tr (Y* - B 


+ Arj + Y^ + A t ^ L D r7 2)A 0 (D a 2fj, L + H ct2 Y l + -D ct 2A 2 l) 


+ 9 [tfl, 2 ^it-R ,2 + 2t^ gAxt/j^ + >, 


(19) 


where A$, A-i,tR ,2 and t R3 are independent of T R and T R . 


Proof First, we observe that (18) can be written as: 


(v^cr) N(I+1) x (V2^o p ) 2N x exp | - -L [n L tr ii L + fi R tr fi R ] 


2 o 2 


N 


YffY L + Y%Y r + £)(Y£ - B n T 0 )* r (Y* - B n T 0 ) 


T%(D a 2 + D a 2 p + —A r )T l - 2(n L tr D a 2 p + Y R D a 2 + A^ L D a 2)T L 


+^2 T R At LR T L + T %{D G 2 + D a% + - 2 (n R tr D a 2 p + Y%D a 2 + A % R D a 2)T R 


/ / “*l 

Lif 

>Tr Jt l 

l 2 L 


To marginalize T R and T#, we assume that they are independent Gaussian ran¬ 
dom vectors. We can then use the following standard result: if X ~ J\T p {ii : £), then 
Al(exp(t* r X)) = exp(t tr /x + |t* r St). 

Therefore, by integrating first with respect to T^, the marginal likelihood of 6 and 
Tfl is proportional to the product of a factor that is independent of Tl and the following 
term 
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■D„ 




exp(t^' 1 T L )dT L , 


where 


t%i = ^ L tr D <rl +Y t £D a , 

1 xN 


A %L 


d**)-- 2 t 


tr a 
R A 


tr 

LR • 


It is now convenient to define A^ -1 := (D a 2 + D a 2 + Aj Ar). The marginal likelihood 
of G and Tr is proportional to the product of a factor that is independent of Tr and 
the term (27r) JV/ ’ 2 |Ao| 1 / 2 exp { } . 

Therefore, the marginal likelihood of G can be explicitly written as 


(V2n<j)- N ^ x (V^a p )~ 2N x (2 7 r) JV / 2 |A 0 | 1/2 


1 


1 


x exp j 2^2 (MlVl + MrVr) - 2^2 


AT 


Y^Yr + Y^Yr + E( Y n - B"T 0 ) tr (Yi - 


^ exp i ^tc/iAotcp - * 


Tr(-D<t 2 + ^ 2 + Ar)T* - 2(fi B tr D a 2 + Y + A% 
p a z p 


The entire last expression is equal to the product of a term that is independent of 
Trj and the following term: 



(Do - 2 + D a 2 + -^Ar) 



2 

ArrAqArr 



expjtflpTiiJdTR, 


where 


i,i = (t*‘R r Dol + Yr-D ct 


A£rA, 




+ Y f r r d t 


■2 + A t 2 r L D r7 2 )A 0 A LR 


If we now define Aj^ 1 := 
respect to Tr, we have 


'Tr 


1. 


exp ^ — ^TrA^" 1 


= (2tt) jv/ 2 |A 1 | 1/2 exp 


= (D a 2 +D a 2 + AjAr) - (ya) 2 A t [ R A 0 A LR and integrate with 

TrJ exp(tR: i TR)dTR = (2 tt) 7V/:2 |Ai| 1/2 exp j^t^AitRp j 
| 2 (*fl,2 + *R,3)Al (tfl,2 + t«,3) | , 


after evaluating the integral. We finally obtain (19). f 
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5 A Bayesian inference for thermal diffusivity 


In this section, we implement our Bayesian approach to infer the thermal diffusivity #, 
an unknown parameter that appears in the heat equation and measures the rapidity of 
the heat propagation through a material (dos Santos et al. 20051. Temperature data 


are available on the basis of cooling experiments. Synthetic data are used to carry out 
the inference. 


Consider the heat equation (one-dimensional diffusion equation for T(x,t )): 


{ d t T - d x (9(x)d x T) =0, x £ (x L ,x R ), 0 < t < t N < oo 
T(0, t) = T L (t), fe[0,Cv] 

T(l,t) =T R (t), te[0,tjv] 

T{x, 0) = g(x), x £ {x L ,x R ). 


( 20 ) 


We want to infer the thermal diffusivity, 6(x), using a Bayesian approach when 
the temperature is measured at / + 1 locations, Xq = x R , X\, x %,..., Xj-i, Xi = x R: 
at each of the N times, ti, t 2 > • ■ •, tzv- Clearly, this problem is a special case of 0 
where Lg = —d x ( 9{x)d x T ) and 9{x) > 0. We can therefore immediately obtain the 
non-normalized posterior distribution of 9 using the marginal likelihood (19). 


The prior distributions for 9 can be specified in different ways. In this section we 
will consider two cases, when 9 is a lognormal random variable and when 9 depends on 
the space variable x and is modeled by means of a lognormal random field. We focus on 
the second case in Subsection |5.3| We start by discussing the case where the thermal 
diffusivity prior is independent of x. 


If we consider a lognormal prior log# ~ AT (v,t), where v G R. and r > 0, then the 
non-normalized posterior distribution of 6 is given by 


Pi/,r(#|Y 1 ,..., Yn) a 


1 


\/ 7 2/k9' 1 


exp 


(log# — ify 
2^2 


p(Y 1 ,...,Y N |#). (21) 


The posterior distribution of # can be approximated by Laplace’s method (Ghosh et al. 
hapter 4]) to obtain 

P»A9 |Yi.Y n ) 


(2006), [Chapter 4]) to obtain a Gaussian posterior 

1 


27t|U(#)| 


exp {-(# - 9) tr H{9)~ 1 (9 - #)} 


where # is the maximum a posteriori probability (MAP) estimate and H(9) is the 
Hessian matrix of the log posterior evaluated at #. 

To assess the behavior of our method, we introduce a synthetic dataset generated 
with constant #. Let us assume, without loss of generality, that the interval time [0,£jv] 
is equal to [0,1], Xl =0 and x R = 1. 
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Dataset A 

In order to generate data, we solve the initial-boundary value problem for the heat 
equation with Robin boundary conditions: 

d x T{x L ,t ) = - (T(x L ,t) - T out ) > *€[0,1], 

d x T{x R ,t) = ~{T out - T(x R ,t)) , tG [0,1], 

and the initial condition T(x, 0) = T 0 , x G (0,1), where 9{x) = 1 x 10 - ' m 2 /s, h is the 
convective heat transfer coefficient, k denotes the thermal conductivity, - = l(l/m), 
T out = 20 °C and T 0 = 100 °C (see Figure 0. 



Figure 1: Exact solution of the initial-boundary value problem for the heat equation, 
0 = lx 10~ 7 m 2 /s. 

A synthetic dataset (hereafter named dataset A) is generated, with a measurement 
standard error noise of ad = 0.56. 

Before presenting the implementation of our novel technique, we will show how the 
Bayesian method works, using the joint likelihood Q, under the very restrictive as¬ 
sumption that the temperature values at the boundaries are exactly known. 


5.1 Example 1 


Suppose that the thermal diffusivity, 6 , is a random variable with a lognormal prior, 
log# ~ N{v,t). In this case, the non-normalized posterior density for 9 is given by 


P„,t(0|Y 1 ,...,Y n )cx 


1 


s/Znd', 


exp 


(log 0 - vf 

2 T 2 


exp 




N 


2<t 2 


n =1 
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where is used here because the boundary data are known exactly. 




Figure 2: Example 1: Comparison between log-likelihoods (on the left) and log- 
posteriors (on the right) for 6 using different numbers of observations, N, and different 
values of er. 


Given that 6 is a lognormal random variable with v = t = 0.1, the resulting posterior 
will depend on a and the number of observations N that are used to compute the log- 
likelihood. Figure [2] shows the behavior of the log-likelihood and the log-posterior for 6 
using different values for a and N. 


80 

60 

40 

20 

0.8 0.9 1 1.1 1.2 1.3 1.4 

Thermal diffusivity 0(1 x 10~ 7 m 2 /s) 

Figure 3: Example 1: Lognormal prior and approximated Gaussian posterior densities 
for 6 1 where a p = a = 0.5 and N = 60. 
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We then use the Laplace approximation to derive the Gaussian posterior approxi¬ 
mated density for 8. The prior and posterior densities for 8 are presented in Figure [3j 
where it can be appreciated that we obtained a Gaussian posterior, with mean 1.0025 
and standard deviation 0.0044, which is concentrated around the true value of the pa¬ 
rameter, 8, despite having a very broad prior. We are now in the position to extend our 
implementation to embrace the case where the temperature values at the boundaries 
are unknown parameters as well. 


5.2 Example 2 


In this example, we consider again 8 as a random variable with a lognormal prior, 
log# ~ A f (y,r). Unlike Example 1, we assume noisy boundary measurements and a 
Gaussian prior distribution for the boundary parameters as in 0 Therefore, the 
non-normalized posterior density for 8 is given by 


p I/ , r (#|Y 1 ,...,Y N ) a 


\/2n8', 


exp 


(log @ vf 

2t 2 


p(Y 1 ,...,Y n |0), 


where p(Yi,..., Yn|#) is the marginal likelihood of 8 defined in Theorem 4.1 


Remark 5.1 Since 8(x) is supposed to be constant, we can obtain equation (|15[ > alter¬ 
natively by solving the heat equation using finite differences (see Appendix C\b 1). 




Figure 4: Example 2: Comparison between log-likelihoods (on the left) and log- 
posteriors (on the right) for 8 using different numbers of observations, N, and different 
values of a. 

Numerical results are now presented using the synthetic dataset A and assuming 
that 8 is a lognormal random variable with v = r = 0.1. Figure [3] shows the behavior of 
the log-likelihood and the log-posterior for 8 using different values of N and cr. Clearly, 
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Figure 5: Example 2: Comparison between log-likelihoods (on the left) and log- 
posteriors (on the right) for 9 using different values of a and a p , with N = 60. 


the accuracy of the estimated 9 depends on the size of the dataset, N , and the reliability 
of measurement devices, cr. 

Figure [5] shows the relationship between the prior distribution of the boundary con¬ 
ditions and their measurements. Although we notice different behaviors of the log- 
likelihood and the log-posterior, these functions exhibit the same argument of the max¬ 
imum which is close to the true value of 9. 


80 
70 
60 
50 
40 
30 
20 
10 

0.8 0.9 1 1.1 1.2 1.3 1.4 

Thermal diffusivity 9(1 x 10~ 7 m 2 /s) 

Figure 6: Example 2: Lognormal prior and approximated Gaussian posterior densities 
for 9 where <j p = a = 0.5 and N = 60. 
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Again, the Laplace approximation is used to derive the Gaussian posterior approx¬ 
imated density for 8. The prior and posterior densities for 8 are presented in Figure [6] 
in which the Gaussian posterior, with mean 0.9955 and standard deviation 0.0047, is 
concentrated around the true value of 8 unlike the very broad prior. 


Information Divergence and Expected Information Gain 


In the Bayesian setting that we adopted to infer the thermal diffusivity, 8 , the utility 
of the performed experiment, given an experimental setup £, can be conveniently mea¬ 


sured by the so-called information divergence (or discrimination information as Kullback 
(1987) called it), which is here defined as the Kullback-Leibler divergence (Kullback and 


Leibler (1951)) between the prior density function p(8) and the posterior density func¬ 
tion Y 1 ,...,Y n ,0: 


D K l( Yu-.Yn.O := / log 


p(fl|Y 1 ,...,Y N ,Q 

p(8) 


p(8 |Y 1 ,...,Y n ,0^. (22) 


The quantity in (22) is always non-negative; it is equal to zero when the prior 


and the posterior coincide; it provides a quantification of the relative discrimination 
between the prior and the posterior; and it depends on the observations Yi,..., Yn. 
Therefore, given the synthetic dataset, A, we may introduce different experimental 
setups of interest by varying the interval time during which the temperature is measured. 
By choosing some specific thermocouples,we may evaluate the information divergence 
for any experimental setup. Moreover, under the same generating process used for the 
dataset A, we may obtain as many synthetic datasets as needed to explore the properties 
of the proposed simulated experiment. The utility of such computer-based experiments 


can be adequately summarized by the so-called expected information gain (Long et al. 


(2013)), which is defined as the marginalization of Dkl over all possible simulated data: 


m ■■= 


log 




p(0|y 1 ,...,y n ,q 

p(8) 


P(8 |Y l5 ..., Y n , Oddp({ Y ; }f =1 |0d({Y ; }f =1 ). 

(23) 

This quantity (231 provides a criterion to determine which features of the setup, £, are, 
on average, most informative when inferring 8. A larger value of /(£) when, say, £ £ A, 
suggests that, given the proposed statistical model, the inference on the unknown pa¬ 
rameter will be more efficient, on average, when the features of the designed experiment 
take value in the set A. 


Let us label as TCI,..., TC7 the thermocouples from the left boundary to the right 
boundary, respectively. 

The numerical estimations of the information divergence for the synthetic dataset 
A and of the expected information gain, by using © to compute the approximated 
posterior in Example 1, are shown in Figures 0® and [9j which refer to the following 
three experimental setups (es’s): 
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esl) £ consists of three non-overlapping time intervals, with the same length, which 
cover the entire observational period [0,1]; 

es2) £ consists of the five inner thermocouples; 

es3) £ is the combination of the two previous experimental setups, esl and es2. 



1 - 1 - 




Information Divergence 
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. 0,,\ 

" Expected Information Gain 
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\ 
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1 


First Time Interval Second Time Interval Third Time Interval 


Figure 7: Example 2: The expected information gain compared with the information 
divergence for the synthetic dataset, A, for the three time intervals experimental setup 
(esl). 

From the values in Figure [TJ which depicts the results from an experimental setup 
in which the temperature measurements are collected at different time intervals, we 
may conclude, by virtue of the interpretation of the expected information gain, that the 
second time interval is the most informative time interval from which to draw inferences 
on the thermal diffusivity, whereas the last time interval is the least informative one. 

Figure [8] summarizes how the expected information gain behaves for the five inner 
thermocouples experimental setup (es2). Given the synthetic dataset A, the sixth ther¬ 
mocouple (TC6) is the one where the information divergence takes the smallest value. 
However, when we look at the expected information gain, we may appreciate the nearly 
symmetric informative content of the thermocouples with respect to the central thermo¬ 
couple (TC4) and how the expected gain about the thermal diffusivity becomes larger 
near the central thermocouple. 

Finally, we look for the best combination of the two previous experimental setups, 
and the corresponding results are displayed in Figure [9] We observe that the highest 
expected information gain is attained at the middle thermocouple (TC4) by using the 
information collected during the second time interval. Any indication provided by the 
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Figure 8: Example 2: The expected information gain compared with the information 
divergence for the synthetic dataset, A, for the five inner thermocouples experimental 
setup (es2). 



Figure 9: Example 2: The expected information gain computed for the combination 
(es3) of the three time intervals and the five inner thermocouples experimental setups. 


numerical estimation of the expected information gain is very valuable to an experi¬ 
mentalist, since it suggests the most relevant features to be considered to build up an 
efficient experiment to infer the unknown parameters of the assumed statistical model. 
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Predictive Posterior Distribution 


In this section, we examine the possibility of predicting the observable temperature 
at future time intervals after estimating the thermal diffusivity, 9. More specifically, 
assume we have inferred 9 using temperature measurements up to time t n . Then, 
we want to predict the temperature in the next time step, t n +\. Given our Bayesian 
model, it is necessary to assume the knowledge of the boundary temperature at time 
t n - )-i. A typical situation could be given by an experiment in which there is interest in 
temperature values at inner points for different boundary values. 

The predictive posterior distribution, p(Y n+ i| {Y k }]( =1 , Tr^+i), is given by 


[ p{ Y n+ 1 |{Y k }£ =1 
Je 

and it is estimated by averaging 


-^~L,n+ 1 1 1 5 0)p{9\ {Y k }™ = i)d 6 », 


1 

M 


M 

^p(Y n+1 |{Y k }" =1 ,r L ,n+l> TR,n+\i 6 i) 
i =1 


where the 9i s are sampled from the posterior distribution of 9. 


(24) 


Figure [TO] shows the one-step-ahead predictive posterior densities at three different 
inner thermocouples based on the observations until time t = 0.5, when the observed 
temperature at thermocouples TC2,TC3 and TC4 were 53.98,55.53 and 57.84°C re¬ 
spectively. 

We emphasize that this methodology allows us to obtain the k-step ahead predictive 
posterior density for any inner thermocouple, assuming boundary conditions subject to 
uncertainty that are adequate for the experiment. 


5.3 Example 3 

In this example, we consider the case where the thermal diffusivity depends on the space 
variable, x. The finite element method used to solve the heat equation under such an 
assumption was presented in Section [3] 


Prior distribution of 9(x) 

Assume that the prior distribution of 9(x) is a lognormal random field with a squared 
exponential (SE) covariance function. Then, the prior distribution of log0(a;) can be 
expressed using the joint multivariate Gaussian distribution: 

(log( 6 »(a"i)),..., log( 6 >(cc s ))) ~ N s (n, I \), (25) 


where fi = (/x,/x,...,/x) tr , K zj = Cov{\og{9(x i )),\og(9(x j ))) = rj 2 exp 1 ) ,i,j = 

77 is the magnitude, and i denotes the length scale. 
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Figure 10: Example 2: The predictive posterior densities of the observable temperatures 
for the thermocouples TC2,TC3,TC4 at time t = 0.52. 


We assume the following priors for the hyperparameters /i, 77 and t (the prior density 
of p and 77 is displayed in Figure [TTj) : 

/x ~ AT (0.1,0.1), 77 ~ half-Cauchy (0.1), £ ~ U (0.5,5). 


In this example, we choose a Gaussian prior for p and uninformative uniform prior 
for l. The half-Cauchy prior for rj was chosen because it is a practical prior for scale 
parameters in hierarchical models (Poison and Scott||2012) (Gelman||2006). 



Figure 11: Example 3: Joint prior density for the hyperparameters p and 77 , with 
p ~ Af (0.1,0.1) and 77 ~ half-Cauchy (0.1). 
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Joint posterior distribution of the hyperparameters p,,r],£ 

Given 6 := (9(x i),..., 0(x s )) tr , let us consider the joint posterior density of the hyper¬ 
parameters (p,,rj,£) that characterize the distribution of logP(:r). 


p(Mi» 7 >^|Yi,...,Yn) oc p(n,r},£) [ p(0\p,r],£)p(Y 1 ,... ,Y N \9)d0. 

•/© 

Let us introduce the auxiliary variable z = (z ±,..., z s ) ~ Af s ^0, C = ^ A' j and consider 

the change of variables transformation: log( 0 j) = p + pzi where Qi := 9(xi ), Zi := z(xi), 
i = 1,..., s. Then, the prior density of 6 is given by 




9\ 62 • • ■ 9 S 
(2 7 r)-i|C , |-5 

r jS e s^+ri(zi+...+z s ) 

p{z\£) 


(log 9 — n) tr C {log 6 - p,) 
2 rf 


exp ( — -z C z 


^SgSH+ritzi+.-.+Zs) ' 

The posterior density of the hyperparameters can be therefore written as 

p(z\£) 


p(p,t,V |Yi,...,Y n ) oc p(p,f,7]) J 


\J\ p(Yi, - - -, Y N |/z, 77, z)dz, 


/ z ? yS e SM+^(zi + ...+z,) 

where J is the Jacobian matrix of the transformation. 

By considering £ as a nuisance parameter, we obtain 

p(Pi f?|Yi,..., Y n ) oc p(p, rf) f p(t) [ p{z\£)p(Y ll ...,Y- N \p 1 rpz)dzd£ (26) 

Ji J z 

after t is marginalized. 

To evaluate the posterior distribution of the hyperparametrs, we need to compute 
the s + 1 dimensional integral in formula (26 1 . Alternatively, Monte Carlo method can 
be used to approximate these integrations. First, we sample i from its prior distri¬ 
bution, p{£). Then, given £ we can sample z and evaluate the joint likelihood func¬ 
tion, p(Yi,..., Yn|/u, 77 , z), for any pair Therefore, we approximate the non- 

normalized posterior distribution using a double sum as follows: 

^pW J p( z K)p(Yi, ..., Y N |/x, r), z)dzd£ « ^ J p(z|^)p(Yi,..., Y N |/z, r],z)d‘. 


i=1 
1 1 

Mi M z 


f z 

Mi M z 


EEp( Y 1 ,.-.,Y N | Ai ,» 7 ) z j ), 


»=1 J=1 


where Zj ~ p(z|^). 
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Figure 12: Example 3: Non-normalized joint posterior density of the hyperparameters 
p and p. The maximum a posteriori probability (MAP) estimate is (—0.05,0.025). 



H T] 


Figure 13: Example 3: Laplace’s approximation of the posterior density of the hyper¬ 
parameters p and p. 


Figure 12 shows that the non-normalized posterior of the hyperparameters p and 
p has a unique mode at (—0.05,0.025). We use then Laplace’s method to obtain a 
Gaussian posterior, using the synthetic dataset A, as shown in Figure [l3| 

A new dataset is now introduced to test our method when 9 depends on x. 

Dataset B 

To analyze the performance of our inferential technique in the case where the thermal 
diffusivity parameter depends on the space variable, x, we consider another synthetic 
dataset (hereafter named dataset B) that is generated similarly to the dataset A, except 


for the fact that 9(x) is sampled randomly from the new prior (25) where p = 0, p = 0.1 
and £ = 5. 


Again, we approximate the posterior distribution of the hyperparameters p and p 
using Laplace’s approximation given the following priors for the hyperparameters p , p 
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and £: 

H ~ J\f (0,0.25), 77 ~ half-Cauchy (0.5), £~t/(4,6), 

where we assume broad priors for 77 and 77 with a more informative unifrom prior for £. 



Figure 14: Example 3: Joint prior density for the hyperparameters 77 and 77 , with 
/i ~ Af (0, 0.25) and 77 ~ half-Cauchy (0.5). 



Figure 15: Example 3: Non-normalized joint posterior density of the hyperparameters 
77 and 77 . The maximum a posteriori probability (MAP) estimate is (0.05,0.025). 

From Figure [l5| we find that the maximum a posteriori probability (MAP) estimate 
is (0.05, 0.025) and Laplace’s approximation can be used. By comparing the prior and 
posterior densities for /i and 77 in Figures |T4| and fl 6 j we can say that the experiment is 
informative since the posterior concentrates around (0.05,0.025) which is close to the 
true value. 


6 Conclusion 

In this work, we developed a general Bayesian approach for one-dimensional linear 
parabolic partial differential equations with noisy boundary conditions. First, we de- 
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H T] 


Figure 16: Example 3: Laplace’s approximation for the posterior density of the hyper¬ 
parameters p and rj. 


rived the joint likelihood of the thermal diffusivity 9 and the boundary parameters. 
Second, we approximated the solution of the forward problem, by showing that such 
solution can be written as a linear function of the boundary conditions. After that, 
we marginalized out the boundary parameters, under the assumptions that they are 
well approximated by piecewise linear functions and that they are independent Gaus¬ 
sian random vectors. This approach can be generalized to any well-posed linear partial 
differential equation. 

On the implementation side, we computed the log-posterior of the thermal diffusiv¬ 
ity in different cases. Besides, we used the Laplace approximation to obtain a Gaussian 
posterior. In the first example, we used directly the joint likelihood of the thermal dif¬ 
fusivity 9 and the boundary parameters, assuming that the boundary conditions were 
known. In the second example, we used the marginalized likelihood of 9, assuming 
that 9 is a lognormal random variable and, as in the previous example, we obtained an 
approximated Gaussian posterior distribution, showing that the unknown value of the 
thermal diffusivity is recovered almost exactly. Moreover we explored two important 
advantages of using the Bayesian approach, by providing the estimation of the expected 
information gain for different experimental setups and the predictive posterior distri¬ 
bution of the temperature. We noticed that the temperature measurements from the 
middle thermocouple at the second time interval are, in general, the most informative 
measurements. Finally, we considered the case where 9 is a lognormal random field 
with squared exponential covariance function. In this case, we obtained the joint poste¬ 
rior distribution for the covariance hyperparameters by applying hierarchical Bayesian 
techniques. 
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Appendix A 

Proof of Theorem 13.21 

First, let us introduce the vectors 
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and the matrices 
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Then, the solution of © is given by: 

u„+i = Bu n + (M + A tSg) 1 ( Pl,uLl + Fji,nTjj) . 

Applying recursively the previous relation we derive the discrete representation 
(Duhamel’s formula) 

n 

U n = B"u 0 + ^ B”- fc (M + A tSe)- 1 (F L , fe T L + F*,*T fl ). 
fe =1 


Now we can build the matrices A n (0),Ai, !n (9) and A R , n (6) ,n = 1,..., N, intro¬ 
duced in the expression (14), to recover the solution of the problem © as a linear 
function of the initial-boundary conditions, namely: 


A n {9) 

= B", 

Al,u{9) 

II 

Cd 


k= 1 

A R ,n(6) 

ii 

td 


L,k, 


and 


Appendix B 

Proof of Theorem 13.31 


From and (14), we can write: 


T n = B' l u 0 + A L ^ n {9)T L + A Rin {9) T R — T L , n F L ^ — T R>n F R p, 
Now, define Al jW ( 0) and A R , n (9) by: 


Al,u{ 9) = A Ltn (9) + 


B n F L - 


0 




0 


(7-l)xl 1) X (rj.—!) (/_ 1) X 1 (/-l)x(JV-n) 


A R , n {0) = A R , n (G)+l B ’ lF «4 


0 


—F-, 


R, 1 


0 


(/-l)xl ( J —!)x (n—!) (;_i)xl (J-l)x(AT- 


Therefore, we obtain equation (15). 


Appendix C 


Proof of Remark 15.1 


Consider the following backward Euler discretization of the local problem ([3]) in the 
interval time, {t n ,t n+ 1 ) = (nAt, (n+ l)At): 


-£i(Ti t n+ 1— 7i >n ) — -^2 (Ti+i t „+i — 2Tj )n _|_i + Tj_i in _)_i) — 0, * — 2, ...,J 

T L ,n = T L (nAt ), 

T R , n = T R {nAt). 


( 27 ) 





F. Ruggeri, Z. Sawlan, M. Scavino and R. Tempone 
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To write the discretization (27) in a matrix form, let us introduce the vectors T„ = 

(i- i)xi 

(T 2 ,„,..., Tj tU ) tr , n = 1,..., N, and the matrix 

/ -2 1 0 0 ... 0 \ 

1-2 1 0 ... 0 

0 1 -2 1 ... 0 

A 

(/—l)x( 7 — 1 ) 


0 

V o 


... 0 

0 0 


1 -2 1 

... 1 -2 


In this way, we may write 


At (T n+ i T n ) Aa , 2 AT„+i Aa , 2 


(Ti,n+lV + Tr^+iw) , 


(28) 


where v = (1, 0,..., 0) tr and w = (0,..., 0, l) tr . 

(/-l)xl (/-l)xl 

The expression (281 is equal to 

(Ii— l ~ 0 a ^-A)T„_|_ i = T n + 9 A ^ 2 |-i v + Tr^+iw) 


and letting = A and B = (/j_i — 9 AA) - 1 , we obtain 

T m +i = BT„ + 9\(TL tTl+ iQ'v + T Rn+ iBw). 


Applying recursively the previous relation, we derive 

n n 

T n = B"T 0 + 9X Y T Ltk B"- fc+1 v + 9\ Y T R , k B n - k+1 w , 

k -1 fc =1 

whose compact matrix form is 

T n = B”T 0 + C„T L + D„Tr , 

where 

• T l = (T Lj1 , .. •, T L n ) tr = (Ti(Ai),... ,T L (7rAf)) tr , 

nx 1 

. Tr = (Trp, ..., TR jn ) tr = (Tfl(At),..., Ti?(nAf)) tr , 

nx 1 

• the matrix C n has column vectors c*. = 0AA" _fc+1 v , fc = 1,..., n , 

(/— 1 ) Xn 

• the matrix D„ has column vectors = 0AA' l-fc+1 w , k = 1,..., n. 

(/— 1 ) Xn 
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Now, we can build the matrices Al u (9) and Ar u { 6) ,n = 1 introduced in 

(I—l)xN (I-l)xN 

the expression (161, to recover the solution of the problem ([3| for each interval time, 
(i n _ i, t n ) , n = 1,..., iV , as a linear function of the initial-boundary conditions: 


A L , n (0) = 

(J-l)xJV 


Cn 

(/— l)xn 


AR, n (9) = 

(I-l)xN 


c n 

(I-l)xn 


0 ") 

(/—l)x (N—n) ) ’ 


0 ") 

(I— 1) X (N—n) ) ■ 


I 



